****Figure 1
use "Figure 1 Data.dta", clear

collapse (sum) slavepop, by(censusyr state_type)
replace slavepop=slavepop/1000
reshape wide slavepop, i(censusyr) j(state_type)

 
twoway (area  slavepop3 censusyr, lpattern(longdash) ) ///
(area slavepop2 censusyr) ///
(area slavepop1 censusyr, ///
legend(label(3 "North") label(2 "Border") label(1 "South") rows(1)) ///
 ytitle("Slave Population" "(1,000's)") xtitle(Census Year) scheme(s2mono) ///
color(black) lpattern(dash) graphregion(color(white)))
 
 
 **figure 2
  
use "Data.dta", clear
gen axis=_n + 0 in 1/19
gen mean_s=.
gen ll_s=.
gen ul_s=.
gen mean_n=.
gen ul_n=.
gen ll_n=.

forvalues c=1(1)19{
summ perslv if south==1 & cong==`c'
replace mean_s=r(mean) if axis==`c'
replace ll_s=r(mean)-0.5*(r(sd)) if axis==`c'
replace ul_s=r(mean)+0.5*(r(sd)) if axis==`c'
summ perslv if south==0 & cong==`c'
replace mean_n=r(mean) if axis==`c'
replace ll_n=r(mean)-0.5*(r(sd)) if axis==`c'
replace ul_n=r(mean)+0.5*(r(sd)) if axis==`c'
}
twoway (rspike ll_n ul_n axis) (scatter mean_n axis) (rspike ll_s ul_s axis) ///
(scatter mean_s axis, scheme(s2mono) legend(order( 2 "Avg. Non-South" 4 "Avg. South" ///
  1 "Standard Deviation") )  xscale(range(1 19)) xtitle("Congress") ///
  xlabel( 1 5 10 15 19) ytitle("District Slave Percentage") graphregion(color(white)))
  
  
***********Figure 3  
use "Data.dta", replace
gen axis=_n+0 in 1/19

gen diff_pct=.
gen p_value_pct=.
forvalues c=1(1)19{
ttest pctvote if cong==`c', by(south)
replace diff_pct=r(mu_2)-r(mu_1) if axis==`c'
replace p_value_pct=r(p) if axis==`c'
}
twoway (line diff_pct axis) (scatter diff_pct axis if p_value_pct<=0.05, mcolor(black)) ///
(scatter diff_pct axis if p_value_pct>0.05, scheme(s2mono) legend(order(1 "Difference" ///
2 "Stat. Sig.") size(small) ) xtitle(Congress, size(small)) xscale(range(1 19)) ///
xlabel( 1 5 10 15 19) ytitle("Pct. Vote", size(small)) yline(0, lcolor(gs8) ///
 lpattern(dash)) name(g1, replace) graphregion(color(white)) ///
)


gen diff_enep=.
gen p_value_enep=.
forvalues c=1(1)19{
ttest enep if cong==`c', by(south)
replace diff_enep=r(mu_2)-r(mu_1) if axis==`c'
replace p_value_enep=r(p) if axis==`c'
}
twoway (line diff_enep axis) (scatter diff_enep axis if p_value_enep<=0.05, mcolor(black)) ///
(scatter diff_enep axis if p_value_enep>0.05,scheme(s2mono) legend(order(1 "Difference" ///
2 "Stat. Significant" 3 "Stat. Insignificant") ) xtitle(Congress, size(small)) xscale(range(1 19)) ///
xlabel( 1 5 10 15 19) ytitle("Eff. N of Parties", size(small)) yline(0, lcolor(gs8) ///
 lpattern(dash)) name(g2, replace) graphregion(color(white)) ///
)

gen diff_enec=.
gen p_value_enec=.
forvalues c=1(1)19{
ttest enec if cong==`c', by(south)
replace diff_enec=r(mu_2)-r(mu_1) if axis==`c'
replace p_value_enec=r(p) if axis==`c'
}
twoway (line diff_enec axis) (scatter diff_enec axis if p_value_enec<=0.05, mcolor(black)) ///
(scatter diff_enec axis if p_value_enec>0.05, scheme(s2mono) legend(order(1 "Difference" ///
2 "Stat. Significant" 3 "Stat. Insignificant") ) xtitle(Congress,size(small)) xscale(range(1 19)) ///
xlabel( 1 5 10 15 19) ytitle("Eff. N of Candidates", size(small)) yline(0, lcolor(gs8) ///
 lpattern(dash)) name(g3, replace) graphregion(color(white)) ///
) 

gen diff_comp=.
gen z_score=.
forvalues c=1(1)19{
prtest competitive if cong==`c', by(south)
replace diff_comp=r(P_2)-r(P_1) if axis==`c'
replace z_score=r(z) if axis==`c'
}
gen p_value_comp=2*(1-normal(abs(z_score)))
twoway (line diff_comp axis) (scatter diff_comp axis if p_value_comp<=0.05, mcolor(black)) ///
(scatter diff_comp axis if p_value_comp>0.05, scheme(s2mono) legend(order(1 "Difference" ///
2 "Stat. Significant" 3 "Stat. Insignificant") ) xtitle(Congress, size(small)) xscale(range(1 19)) ///
xlabel( 1 5 10 15 19) ytitle("Prop. of Competitive Districts", size(small)) yline(0, lcolor(gs8) ///
 lpattern(dash))name(g4, replace) graphregion(color(white)) ///
)

grc1leg g1 g2 g3 g4, graphregion(color(white)) legendfrom(g1)

*******************************************************************************
*** Generate Predicted Values for Figure 3 and Estimates for Models 1-4
set seed 12345
*model 1
use "Data.dta", replace

egen district=group(statecode cd)
egen new_state=group(statecode)
recode2 totpop -9999=.
tsset district cong
replace totpop=totpop/1000
mixed pctvote  perslv m totpop   open  /// 6
i.cong || new_state: || district: , residuals(ar 1, t(cong)) 


local a=1
foreach v of varlist perslv m totpop  {
quietly summarize `v' if e(sample)
local val`a'=r(mean)
local sd`a'=r(sd)
local a=`a'+1
}
*
matrix v=e(V)
matrix b=e(b)

tempfile sim1
drawnorm MG_b1-MG_b28, n(10000) means(b) cov(v) clear
postutil clear
postfile mypost model exv ul_95 ll_95 ul_90 ll_90 using `sim1', replace
scalar h_x1_l=`val1'-`sd1'
scalar h_x1_h=`val1'+ `sd1'
scalar h_x2=`val2'
scalar h_x3=`val3'
scalar h_x4=0
* cong 2/
gen x_betahat_l=(MG_b1*h_x1_l + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

gen x_betahat_h=(MG_b1*h_x1_h + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

generate diff_expv=x_betahat_h-x_betahat_l	
_pctile diff_expv, p(2.5, 97.5)
local ll_95 =r(r1)
local ul_95 =r(r2)
_pctile diff_expv, p(5, 95)
local ll_90 =r(r1)
local ul_90 =r(r2)
summarize diff_expv
local exv=r(mean)
local model=1
di `exv'
di `ll'
di `ul'
post mypost (`model')  (`exv') (`ul_95') (`ll_95') (`ul_90') (`ll_90')		
postclose mypost
use `sim1', clear

*model 2
use "Data.dta", replace
egen district=group(statecode cd)
egen new_state=group(statecode)
recode2 totpop -9999=.
tsset district cong
replace totpop=totpop/1000

mixed enep  perslv m totpop open /// 6
i.cong  || new_state: || district:, residuals(ar 1, t(cong))


local a=1
foreach v of varlist perslv m totpop  {
quietly summarize `v' if e(sample)
local val`a'=r(mean)
local sd`a'=r(sd)
local a=`a'+1
}
*
matrix v=e(V)
matrix b=e(b)

tempfile sim2
drawnorm MG_b1-MG_b28, n(10000) means(b) cov(v) clear
postutil clear
postfile mypost model exv ul_95 ll_95 ul_90 ll_90 using `sim2', replace
scalar h_x1_l=`val1'-`sd1'
scalar h_x1_h=`val1'+ `sd1'
scalar h_x2=`val2'
scalar h_x3=`val3'
scalar h_x4=0
* cong 2/
gen x_betahat_l=(MG_b1*h_x1_l + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

gen x_betahat_h=(MG_b1*h_x1_h + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

generate diff_expv=x_betahat_h-x_betahat_l	
_pctile diff_expv, p(2.5, 97.5)
local ll_95 =r(r1)
local ul_95 =r(r2)
_pctile diff_expv, p(5, 95)
local ll_90 =r(r1)
local ul_90 =r(r2)
summarize diff_expv
local exv=r(mean)
local model=2
di `exv'
di `ll'
di `ul'
post mypost (`model')  (`exv') (`ul_95') (`ll_95') (`ul_90') (`ll_90')	
postclose mypost
use `sim2', clear

* model 3
use "Data.dta", replace
egen district=group(statecode cd)
egen new_state=group(statecode)
recode2 totpop -9999=.
tsset district cong
replace totpop=totpop/1000
mixed enec  perslv m totpop open   /// 6
i.cong  || new_state: || district:, residuals(ar 1, t(cong))


local a=1
foreach v of varlist perslv m totpop  {
quietly summarize `v' if e(sample)
local val`a'=r(mean)
local sd`a'=r(sd)
local a=`a'+1
}
*
matrix v=e(V)
matrix b=e(b)

tempfile sim3
drawnorm MG_b1-MG_b28, n(10000) means(b) cov(v) clear
postutil clear
postfile mypost model exv ul_95 ll_95 ul_90 ll_90 using `sim3', replace
scalar h_x1_l=`val1'-`sd1'
scalar h_x1_h=`val1'+ `sd1'
scalar h_x2=`val2'
scalar h_x3=`val3'
scalar h_x4=0
* cong 2/
gen x_betahat_l=(MG_b1*h_x1_l + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

gen x_betahat_h=(MG_b1*h_x1_h + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

generate diff_expv=x_betahat_h-x_betahat_l	
_pctile diff_expv, p(2.5, 97.5)
local ll_95 =r(r1)
local ul_95 =r(r2)
_pctile diff_expv, p(5, 95)
local ll_90 =r(r1)
local ul_90 =r(r2)
summarize diff_expv
local exv=r(mean)
local model=3
di `exv'
di `ll'
di `ul'
post mypost (`model')  (`exv') (`ul_95') (`ll_95') (`ul_90') (`ll_90')		
postclose mypost
use `sim3', clear

***************** model 4

use "Data.dta", replace

egen district=group(statecode cd)
egen new_state=group(statecode)
recode2 totpop -9999=.
tsset district cong
replace totpop=totpop/1000
btscs competitive cong district, gen(time) 

gen time_SQ=time^2
gen time_C=time^3

melogit competitive perslv m totpop   time time_SQ time_C ///
 open i.cong || new_state: || district:

local a=1
foreach v of varlist perslv m totpop  time time_SQ time_C {
quietly summarize `v' if e(sample)
local val`a'=r(mean)
local sd`a'=r(sd)
local a=`a'+1
}
*
matrix v=e(V)
matrix b=e(b)

tempfile sim4
drawnorm MG_b1-MG_b29, n(10000) means(b) cov(v) clear
postutil clear
postfile mypost model exv ul_95 ll_95 ul_90 ll_90 using `sim4', replace
scalar h_x1_l=`val1'-`sd1'
scalar h_x1_h=`val1'+ `sd1'
scalar h_x2=`val2'
scalar h_x3=`val3'
scalar h_x4=`val4'
scalar h_x5=`val5'
scalar h_x6=`val6'
scalar h_x7=0

* cong 2/ is coeff c9
gen x_betahat_l=(exp(MG_b1*h_x1_l + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4+ MG_b5*h_x5 ///
+ MG_b6*h_x6 + MG_b7*h_x7 + MG_b9 + MG_b27))/(1+exp(MG_b1*h_x1_l + ///
MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4+ MG_b5*h_x5 ///
+ MG_b6*h_x6 + MG_b7*h_x7 + MG_b9+ MG_b27))

gen x_betahat_h=(exp(MG_b1*h_x1_h + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4+ MG_b5*h_x5 ///
+ MG_b6*h_x6 + MG_b7*h_x7 + MG_b9 + MG_b27))/(1+exp(MG_b1*h_x1_h + ///
MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4+ MG_b5*h_x5 ///
+ MG_b6*h_x6 + MG_b7*h_x7 + MG_b9 + MG_b27))

generate diff_expv=x_betahat_h-x_betahat_l	
_pctile diff_expv, p(2.5, 97.5)
local ll_95 =r(r1)
local ul_95 =r(r2)
_pctile diff_expv, p(5, 95)
local ll_90 =r(r1)
local ul_90 =r(r2)
summarize diff_expv
local exv=r(mean)
local model=4
di `exv'
di `ll'
di `ul'
post mypost (`model')  (`exv') (`ul_95') (`ll_95') (`ul_90') (`ll_90')		
postclose mypost



use `sim1', clear
append using `sim2'
append using `sim3'
append using `sim4'

format %9.3f exv 
twoway (rspike ul_95 ll_95 model if model<4, yaxis(1) lcolor(gs10)) ///
(rcap ul_90 ll_90 model if model<4, yaxis(1) lcolor(gs6)) ///
(scatter exv model if model<4 , mcolor(black)  msymbol(O) mlabel(exv) yaxis(1)) ///
(rspike ul_95 ll_95 model if model==4,  yaxis(2) lcolor(gs10)) ///
(rcap ul_90 ll_90 model if model==4,  yaxis(2) lcolor(gs6)) ///
(scatter exv model if model==4 , mcolor(black)  msymbol(O) mlabel(exv) yaxis(2) ///
 xscale(range(0.75 4.25))  scheme(s2mono) graphregion(color(white)) ///
 xlabel(1 `" "Model 1" "DV: Pct. Vote" " "' 2 `" "Model 2" "DV: ENEP" ""' ///
 3 `" "Model 3" "DV: ENEC" " "' 4 `" "Model 4" "DV: Competitive" ""', labsize(small)) ///
xtitle(" ") legend(order(3 "Estimated Change" 1 "95% C.I." 2 "90% C.I") row(1) size(small) ) ///
 ytitle("Pred. Probability", size(small) axis(2)) ytitle("Percent", size(small) axis(1)) ) 




********************************************************************************
*Generate Predicted Values and Estimates for Models 5-8 and Figure 5
*******************************************************************************

*model 5
set seed 12345
use "Data.dta", replace

egen district=group(statecode cd)
egen new_state=group(statecode)
recode2 totpop -9999=.
tsset district cong
replace totpop=totpop/1000
mixed pctvote  perslv m totpop   open  /// 6
i.cong if south==1 || new_state: || district: , residuals(ar 1, t(cong)) 


local a=1
foreach v of varlist perslv m totpop  {
quietly summarize `v' if e(sample)
local val`a'=r(mean)
local sd`a'=r(sd)
local a=`a'+1
}
*
matrix v=e(V)
matrix b=e(b)

tempfile sim1
drawnorm MG_b1-MG_b28, n(10000) means(b) cov(v) clear
postutil clear
postfile mypost model exv ul_95 ll_95 ul_90 ll_90 using `sim1', replace
scalar h_x1_l=`val1'-`sd1'
scalar h_x1_h=`val1'+ `sd1'
scalar h_x2=`val2'
scalar h_x3=`val3'
scalar h_x4=0
* cong 2/
gen x_betahat_l=(MG_b1*h_x1_l + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

gen x_betahat_h=(MG_b1*h_x1_h + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

generate diff_expv=x_betahat_h-x_betahat_l	
_pctile diff_expv, p(2.5, 97.5)
local ll_95 =r(r1)
local ul_95 =r(r2)
_pctile diff_expv, p(5, 95)
local ll_90 =r(r1)
local ul_90 =r(r2)
summarize diff_expv
local exv=r(mean)
local model=1
di `exv'
di `ll'
di `ul'
post mypost (`model')  (`exv') (`ul_95') (`ll_95') (`ul_90') (`ll_90')		
postclose mypost
use `sim1', clear

*model 6
use "Data.dta", replace
egen district=group(statecode cd)
egen new_state=group(statecode)
recode2 totpop -9999=.
tsset district cong
replace totpop=totpop/1000

mixed enep  perslv m totpop open /// 6
i.cong if south==1 || new_state: || district:, residuals(ar 1, t(cong))


local a=1
foreach v of varlist perslv m totpop  {
quietly summarize `v' if e(sample)
local val`a'=r(mean)
local sd`a'=r(sd)
local a=`a'+1
}
*
matrix v=e(V)
matrix b=e(b)

tempfile sim2
drawnorm MG_b1-MG_b28, n(10000) means(b) cov(v) clear
postutil clear
postfile mypost model exv ul_95 ll_95 ul_90 ll_90 using `sim2', replace
scalar h_x1_l=`val1'-`sd1'
scalar h_x1_h=`val1'+ `sd1'
scalar h_x2=`val2'
scalar h_x3=`val3'
scalar h_x4=0
* cong 2/
gen x_betahat_l=(MG_b1*h_x1_l + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

gen x_betahat_h=(MG_b1*h_x1_h + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

generate diff_expv=x_betahat_h-x_betahat_l	
_pctile diff_expv, p(2.5, 97.5)
local ll_95 =r(r1)
local ul_95 =r(r2)
_pctile diff_expv, p(5, 95)
local ll_90 =r(r1)
local ul_90 =r(r2)
summarize diff_expv
local exv=r(mean)
local model=2
di `exv'
di `ll'
di `ul'
post mypost (`model')  (`exv') (`ul_95') (`ll_95') (`ul_90') (`ll_90')	
postclose mypost
use `sim2', clear

* model 7
use "Data.dta", replace
egen district=group(statecode cd)
egen new_state=group(statecode)
recode2 totpop -9999=.
tsset district cong
replace totpop=totpop/1000
mixed enec  perslv m totpop open   /// 6
i.cong  if south==1 || new_state: || district:, residuals(ar 1, t(cong))


local a=1
foreach v of varlist perslv m totpop  {
quietly summarize `v' if e(sample)
local val`a'=r(mean)
local sd`a'=r(sd)
local a=`a'+1
}
*
matrix v=e(V)
matrix b=e(b)

tempfile sim3
drawnorm MG_b1-MG_b28, n(10000) means(b) cov(v) clear
postutil clear
postfile mypost model exv ul_95 ll_95 ul_90 ll_90 using `sim3', replace
scalar h_x1_l=`val1'-`sd1'
scalar h_x1_h=`val1'+ `sd1'
scalar h_x2=`val2'
scalar h_x3=`val3'
scalar h_x4=0
* cong 2/
gen x_betahat_l=(MG_b1*h_x1_l + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

gen x_betahat_h=(MG_b1*h_x1_h + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4 ///
+ MG_b6 + MG_b24)

generate diff_expv=x_betahat_h-x_betahat_l	
_pctile diff_expv, p(2.5, 97.5)
local ll_95 =r(r1)
local ul_95 =r(r2)
_pctile diff_expv, p(5, 95)
local ll_90 =r(r1)
local ul_90 =r(r2)
summarize diff_expv
local exv=r(mean)
local model=3
di `exv'
di `ll'
di `ul'
post mypost (`model')  (`exv') (`ul_95') (`ll_95') (`ul_90') (`ll_90')		
postclose mypost
use `sim3', clear

**model 8

use "Data.dta", replace

egen district=group(statecode cd)
egen new_state=group(statecode)
recode2 totpop -9999=.
tsset district cong
replace totpop=totpop/1000
btscs competitive cong district, gen(time) 

gen time_SQ=time^2
gen time_C=time^3

melogit competitive perslv m totpop   time time_SQ time_C ///
 open i.cong if south==1 || new_state: || district:

local a=1
foreach v of varlist perslv m totpop  time time_SQ time_C {
quietly summarize `v' if e(sample)
local val`a'=r(mean)
local sd`a'=r(sd)
local a=`a'+1
}
*
matrix v=e(V)
matrix b=e(b)

tempfile sim4
drawnorm MG_b1-MG_b29, n(10000) means(b) cov(v) clear
postutil clear
postfile mypost model exv ul_95 ll_95 ul_90 ll_90 using `sim4', replace
scalar h_x1_l=`val1'-`sd1'
scalar h_x1_h=`val1'+ `sd1'
scalar h_x2=`val2'
scalar h_x3=`val3'
scalar h_x4=`val4'
scalar h_x5=`val5'
scalar h_x6=`val6'
scalar h_x7=0

* cong 2/ is coeff c9
gen x_betahat_l=(exp(MG_b1*h_x1_l + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4+ MG_b5*h_x5 ///
+ MG_b6*h_x6 + MG_b7*h_x7 + MG_b9 + MG_b27))/(1+exp(MG_b1*h_x1_l + ///
MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4+ MG_b5*h_x5 ///
+ MG_b6*h_x6 + MG_b7*h_x7 + MG_b9+ MG_b27))

gen x_betahat_h=(exp(MG_b1*h_x1_h + MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4+ MG_b5*h_x5 ///
+ MG_b6*h_x6 + MG_b7*h_x7 + MG_b9 + MG_b27))/(1+exp(MG_b1*h_x1_h + ///
MG_b2*h_x2 + MG_b3*h_x3 + MG_b4*h_x4+ MG_b5*h_x5 ///
+ MG_b6*h_x6 + MG_b7*h_x7 + MG_b9 + MG_b27))

generate diff_expv=x_betahat_h-x_betahat_l	
_pctile diff_expv, p(2.5, 97.5)
local ll_95 =r(r1)
local ul_95 =r(r2)
_pctile diff_expv, p(5, 95)
local ll_90 =r(r1)
local ul_90 =r(r2)
summarize diff_expv
local exv=r(mean)
local model=4
di `exv'
di `ll'
di `ul'
post mypost (`model')  (`exv') (`ul_95') (`ll_95') (`ul_90') (`ll_90')		
postclose mypost



use `sim1', clear
append using `sim2'
append using `sim3'
append using `sim4'

format %9.3f exv 
twoway (rspike ul_95 ll_95 model if model<4, yaxis(1) lcolor(gs10)) ///
(rcap ul_90 ll_90 model if model<4, yaxis(1) lcolor(gs6)) ///
(scatter exv model if model<4 , mcolor(black)  msymbol(O) mlabel(exv) yaxis(1)) ///
(rspike ul_95 ll_95 model if model==4,  yaxis(2) lcolor(gs10)) ///
(rcap ul_90 ll_90 model if model==4,  yaxis(2) lcolor(gs6)) ///
(scatter exv model if model==4 , mcolor(black)  msymbol(O) mlabel(exv) yaxis(2) ///
 xscale(range(0.75 4.25))  scheme(s2mono) graphregion(color(white)) ///
 xlabel(1 `" "Model 5" "DV: Pct. Vote" " "' 2 `" "Model 6" "DV: ENEP" ""' ///
 3 `" "Model 7" "DV: ENEC" " "' 4 `" "Model 8" "DV: Competitive" ""', labsize(small)) ///
xtitle(" ") legend(order(3 "Estimated Change" 1 "95% C.I." 2 "90% C.I") row(1) size(small) ) ///
 ytitle("Pred. Probability", size(small) axis(2)) ytitle("Percent", size(small) axis(1)) ) 

 
 

